import numpy as np

def schmidt_spike_removal(original_signal, fs):
    """
    Schmidt 尖峰移除
    参数:
      original_signal: 输入信号
      fs: 采样频率
    输出:
      despiked_signal: 移除尖峰后的信号
      zeroed_time_MAA: 被置零样本的总时间（秒）
    """

    # 确保信号是列向量
    original_signal = original_signal.flatten()

    # 确定窗口大小（500 毫秒）
    windowsize = round(fs / 2)

    # 找到不完整窗口中的样本数
    trailingsamples = len(original_signal) % windowsize

    # 将信号重构为多个窗口
    sampleframes = original_signal[:-trailingsamples].reshape(-1, windowsize)

    # 找到最大绝对振幅 (MAA)
    MAAs = np.max(np.abs(sampleframes), axis=1)

    # 初始化置零时间计数器
    zeroed_samples_count = 0

    # 当有样本大于中值MAA的3倍时移除尖峰
    while np.any(MAAs > np.median(MAAs) * 3):
        # 找到具有最大MAA的窗口
        window_num = np.argmax(MAAs)

        # 找到窗口中尖峰的位置
        spike_position = np.argmax(np.abs(sampleframes[window_num, :]))

        # 找到零交叉点（信号变化符号的点）
        zero_crossings = np.abs(np.diff(np.sign(sampleframes[window_num, :]))) > 1
        zero_crossings = np.append(zero_crossings, 0)

        # 找到尖峰的起点
        if np.any(zero_crossings[:spike_position]):
            spike_start = max(1, np.where(zero_crossings[:spike_position])[0][-1])
        else:
            spike_start = 0

        # 找到尖峰的终点
        zero_crossings[:spike_position] = 0
        if np.any(zero_crossings):
            spike_end = min(np.where(zero_crossings)[0][0], windowsize - 1)
        else:
            spike_end = windowsize - 1

        # 将尖峰置零
        sampleframes[window_num, spike_start:spike_end+1] = 0.0001

        # 记录置零的样本数
        zeroed_samples_count += spike_end - spike_start + 1

        # 重新计算MAA
        MAAs = np.max(np.abs(sampleframes), axis=1)

    # 计算被置零样本的总时间（秒）
    zeroed_time_MAA = zeroed_samples_count / fs

    # 将移除尖峰后的信号重构为单个向量
    despiked_signal = sampleframes.flatten()

    # 将拖尾样本加回到信号中
    despiked_signal = np.concatenate((despiked_signal, original_signal[len(despiked_signal):]))

    return despiked_signal, zeroed_time_MAA
